
# Using the Julia script and data files

This is an explanation of how to use the Julia script and data files available
in the Harvard Dataverse to generate and verify the bounds for Witsenhausen's
problem in [^1].

[^1]: A.J.F. Bekker, O. Kuryatnikova, F.M. de Oliveira Filho, and J.C. Vera,
      Optimization hierarchies for distance-avoiding sets in compact spaces,
      arXiv:2304.05429, 2023, 34pp.

If you have questions or remarks, please contact:

  * Fernando Mário de Oliveira Filho <F.M.deOliveiraFilho@tudelft.nl>
    

## The files in the repository

The files in the repository are contained in a tar archive.  Expanding the
archive will create a folder called `Witsenhausen`.

You can use the Julia module `Witsenhausen` to generate and solve the
semidefinite programs that give upper bounds for Witsenhausen's problem (by
using the `three_point_bound` function), but to check the bounds given in [^1]
you can use the solution files provided in the `data` subfolder.

These solution files always follow the same naming convention, namely:

  * `pn-d-0.dat`: data for the bound in dimension ``n`` using the cone ``Q_1^d``
    and no BQP inequalities.  For example, `p3-6-0.dat` contains the solution
    data for the problem in dimension 3 with ``d = 6``.
    
  * `pn-d-1.dat`: data for the bound in dimension ``n`` using the cone ``Q_1^d``
    and BQP inequalities.  The inequalities are given together with the solution
    data.
    
There are also four *sample files* in the `data` subfolder: `sample-512-6.dat`,
`sample-512-10.dat`, `sample-512-14.dat`, and `sample-512-18.dat`.  These
contain the basis of the invariant ring and the sample used to implement the SOS
constraint for ``d = 6, 10, 14, 18``.


## Using the script

To use the Julia script, enter the `Witsenhausen` folder and start Julia.  Then
enter the Pkg REPL of Julia by pressing `]` and type:

    (@v1.8) pkg> activate .
    (Witsenhausen) pkg> instantiate

Go back to the Julia REPL by pressing backspace, then type:

    julia> Using Witsenhausen
    julia> Using Serialization
    
The `Serialization` package is used to load the solution files.


## Loading and inspecting a solution file

Here is how to load a solution file and get basic problem information:

    julia> data = deserialize("data/p3-10-0.dat");
    julia> problem_info(data)
    Problem information (objective value is not rigorous):
             n (dimension): 3
          d (Q_1^d degree): 10
         Inner product t_0: 0.0
        # BQP inequalities: 0
                 Objective: 0.309297759031180254...

The objective value above is a *tentative* upper bound for ``α₃`` given by
our problem with ``d = 10`` and no BQP inequalities.  It is a tentative bound
since the solution still has to be verified and perhaps slightly changed, as
explained in [^1], so we can be sure to have an upper bound.


## Checking the SOS constraint

As explained in [^1], the SOS constraint is not strictly satisfied by the
solution, but its violation should be very small.  You can get the minimum
eigenvalue of any of the ``Q`` matrices in the constraint (see (25) in [^1]) and
the maximum violation of the SOS constraint by using the `compute_sos_violation`
function:

    julia> data = deserialize("data/p3-10-0.dat");
    julia> compute_sos_violation(data, "sample-512-10.dat")
    [...]
    (7.4907398096...e-09, [1.067496757320885e-56 +/- 3.77e-72])

The function returns the minimum eigenvalue of any of the ``Q`` matrices, in
this case ``7.49… × 10^{-9}``, and the maximum violation of the constraint,
which is ``1.067… × 10^{-56}``.  The minimum eigenvalue is much larger than the
maximum violation, and so we know that the solution can be fixed so that the SOS
constraint is then satisfied.

Here it is important to provide a sample file, as above.  Since the solution in
`data` has ``d = 10``, we use the sample file ``sample-512-10.dat``.  The sample
file is important because the basis of the space of invariant polynomials used
to compute the bound depends on the sample, and so to check the solution we want
to use exactly the same sample that we used to generate the problem.


## Checking the linear constraints

In [^1] it is described how to check that a solution indeed satisfies all linear
constraints and, if this is not the case, how to change it so that it does.
This is done by the `fix_linear_constraints` function.  It verifies the linear
constraints, fixes the solution if needed, and shows the objective values before
and after the fix:

    julia> data = deserialize("data/p3-18-1.dat")
    julia> fix_linear_constraints(data);
    [...]
    Result (round last digit up after truncation):
        Original objective: 0.297741420070076956... (rigorous upper bound)

The final objective value given is a rigorous upper bound; if you truncate it,
then you have to round it up.  The function also returns this upper bound.

In this case, all constraints were satisfied by the original solution to the
desired tolerance, so the original objective value is a rigorous upper bound.
Sometimes this is not the case.

If a solution uses many BQP inequalities, then this verification procedure can
take a long time — in the order of minutes, not hours.
